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(2)  Objectives:  Briefly  summarize  the  objectives  of  the  research  effort  or  the 
statement  of  work. 

Statement  of  problem 

We  conjecture  that  a  deterministic  or  random  waveform  that  is  sampled  at  a  rate  less  than 
the  classical  Nyquist  rate  may  be  successfully  reconstructed.  Intuitively,  additional 
information  must  be  available  at  the  sampling  instants,  in  order  to  remove  the  aliased 
spectral  components.  For  example,  the  additional  information  may  be  one  or  more  higher 
order  derivatives.  Another  possibility  is  a  finite  pulse  that  retains  signal  information 
during  the  entire  pulse  duration  samples  the  signal.  Yet  another  possibility  is  that  two  or 
more  closely  spaced  samples  are  retained  each  sampling  instant. 

Our  hypothesis  is  that  a  signal  maybe  sampled  at  one-half  the  Nyquist  rate  if  two 
arbitrarily  closely  spaced  samples  are  retained  each  sampling  instant.  In  this  case  the 
system  sampling  rate  could  be  cut  in  half,  at  the  price  of  carrying  two  samples  per 
sampling  instant.  There  appears  to  be  no  theoretical  reasons  that  the  sampling  rate  cannot 
be  further  reduced. 

Impact  of  the  results  of  research 

A  primary  goal  in  radar  waveform  design  is  to  achieve  a  prescribed  discrete  time  radar 
ambiguity  function.  Radar  waveform  design  to  meet  a  given  ambiguity  function 
specification  is  an  important  issue  under  the  waveform  diversity  technology.  Much  of  the 
development  in  this  area  is  for  continuous  time  signals.  For  the  discrete  time  case,  the 
sampling  rate  becomes  an  important  design  issue.  On  the  one  hand,  too  low  of  a  sampling 
rate  results  in  spectral  aliasing,  on  the  other  hand,  a  sampling  rate  chosen  higher  than 
necessary  increases  the  computational  burden.  We  show  in  this  research  project  that 
aliased  spectra,  arising  from  sampling  below  the  Nyquist  rate,  may  be  completely 
eliminated  in  the  absence  of  quantization  noise  and  timing  errors  (jitter  between  doublet 
samples). 

For  the  purpose  of  this  research,  we  assume  that  the  desired  radar  waveform  is  available. 
Since  radar  processing  is  done  digitally,  it  is  important  to  examine  the  effects  of  sampling 
the  ambiguity  function  in  the  Delay-Doppler  plane  and  study  the  resulting  resolution 
trade-offs,  reconstruction  issues,  and  aliasing  problems.  The  preliminary  techniques 
developed  in  this  research,  and  applied  to  the  ambiguity  function  in  this  proposed 
research,  should  be  quite  useful  in  this  context. 


Approach 


The  approach  for  this  research  project  is  to  apply  the  proposed  aliased  signal 
reconstruction  algorithms  to  radar-like  waveforms  (i.e.  random)  then  form  the  ambiguity 
function  from  the  reconstructed  waveform.  We  will  show  that,  based  on  simulation 
results  and  mathematical  derivation,  that  the  reconstruction  algorithm  is  valid  for 
bandlimited  random  processes. 

In  addition  to  analytically  proving  convergence  of  the  reconstruction  algorithm,  we 
present  plots  using  Matlab-based  simulation  tools  to  illustrate  the  reconstruction  errors 
when  quantization  noise  and  timing  jitter  are  present.  Both  quantization  noise  and  timing 
jitter  will  be  real-world  problems  a  radar  design  engineer  would  need  to  consider  if  these 
algorithms  are  used  in  on-line  processing.  For  off-line  waveform  design  purposes,  noise 
sensitivities  may  be  less  of  an  issue. 


(3)  Status  of  effort:  We  have  achieved  our  planned  objectives  to  derive 
reconstruction  equations  to  remove  the  aliasing  when  sampling  (doublet  impulse 
sampling)  at  Vi  the  Nyquist  rate  (Appendix  A).  The  reconstruction  equations  have 
been  demonstrated  as  valid  for  a  random  waveform  by  both  simulation  and 
mathematical  proof  (Appendix  D).  The  effects  of  quantization  noise  and  jitter  have 
been  illustrated  subjectively  over  several  dozens  test  cases  (Appendix  A).  The 
Matlab  M-files  developed  may  be  found  in  Appendix  C.  An  interim  paper 
presented  at  The  Defense  Applications  of  Signal  Processing  Workshop, 
www.dasp.ws,  may  be  viewed  in  Appendix  B. 

During  the  course  of  this  research  project  we  identified  a  promising  avenue  to 
mitigate  jitter  noise.  This  research  has  not  been  pursued  as  part  of  this  project, 
however,  the  basic  idea  is  documented  in  the  Appendix  A  conclusions  section. 

(4)  Abstract: 

A  primary  goal  in  radar  waveform  design  is  to  achieve  a  prescribed  discrete  time  radar 
ambiguity  function.  Radar  waveform  design  to  meet  a  given  ambiguity  function 
specification  is  an  important  issue  under  the  waveform  diversity  technology.  Much  of  the 
development  in  this  area  is  for  continuous  time  signals.  For  the  discrete  time  case,  the 
sampling  rate  becomes  an  important  design  issue.  On  the  one  hand,  too  low  of  a  sampling 
rate  results  in  spectral  aliasing,  on  the  other  hand,  a  sampling  rate  chosen  higher  than 
necessary  increases  the  computational  burden.  We  show  in  this  project  that  aliased 
spectra,  arising  from  sampling  below  the  Nyquist  rate,  may  be  completely  eliminated.  We 
plot  the  effects  of  quantization  noise  and  timing  jitter  for  subjective  evaluation. 
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Appendix  A:  Technical  Summary  of  Research  Project  “Radar  Waveform  Design  of 

Undersampled  Signals” 


Background 

A  primary  goal  in  radar  waveform  design  is  to  achieve  a  prescribed  discrete  time  radar 
ambiguity  function.  Radar  waveform  design  to  meet  a  given  ambiguity  function 
specification  is  an  important  issue  under  the  waveform  diversity  technology.  Much  of  the 
development  in  this  area  is  for  continuous  time  signals.  For  the  discrete  time  case,  the 
sampling  rate  becomes  an  important  design  issue.  On  the  one  hand,  too  low  of  a  sampling 
rate  results  in  spectral  aliasing,  on  the  other  hand,  a  sampling  rate  chosen  higher  than 
necessary  increases  the  computational  burden.  We  show  in  this  research  project  that 
aliased  spectra,  arising  from  sampling  below  the  Nyquist  rate,  may  be  completely 
eliminated. 

The  radar  ambiguity  function  based  waveform  design  is  complicated  by  the  fact  that  for  a 
given  waveform,  the  ambiguity  function  can  be  readily  calculated.  Flowever,  given  an 
ambiguity  function  specification,  it  is  possible  to  have  more  than  one  waveform  that 
meets  the  specification.  For  the  purpose  of  this  research  project,  we  assume  that  the 
desired  waveform  is  available.  Since  radar  processing  is  done  digitally,  it  is  important  to 
examine  the  effects  of  sampling  the  ambiguity  function  in  the  Delay-Doppler  plane  and 
study  the  resulting  resolution  trade-offs,  reconstruction  issues,  and  aliasing  problems.  The 
preliminary  techniques  developed  in  this  project,  and  applied  to  the  ambiguity  function  in 
this  proposed  research,  should  be  quite  useful  in  this  context. 

The  goal  of  this  research  project  was  to  apply  the  preliminary  aliased  signal 
reconstruction  to  radar  waveform  design  against  ambiguity  plane  constraints  [19].  The 
radar  receiver  discrete  time  matched  filter  computational  complexity  may  be  potential 
reduced  by  implementing  the  signal  restoration  algorithm  summarised  in  this  report.  We 
originally  hypothesised,  based  on  early  simulation  results,  that  the  reconstruction 
algorithm  is  valid  for  bandlimited  random  processes.  Subsequently,  we  proved 
analytically  that  exact  reconstruction  of  bandlimited  random  processes  is  possible.  This  is 
an  important  result  since  modern  radar  waveform  design  involves  the  use  of 
pseudorandom  coding  for  pulse  compression  and  LPI/LPD  waveforms. 

Sampling  Theory 

A  generalised  version  of  the  Nyquist  sampling  theorem  admits  sampling  at  an  average 
rate  equal  to  twice  the  highest  frequency  component  of  the  sampled  signal.  For  example, 
we  may  envision  a  sampling  structure  wherein  a  pair  of  closely  spaced  impulses  perform 
the  signal  sampling,  each  impulse  pair  sampling  the  signal  at  one-half  the  minimum  rate 
required  for  ideal  impulse  sampling.  In  this  report  we  derive  a  frequency  domain 
restoration  algorithm  and  its  corresponding  time  domain  interpolation  formula,  and 
analytically  quantify  the  effect  of  impulse  pair  spacing  on  the  numerical  conditioning  of 


the  restoration  algorithm.  Simulations  results  are  provided  to  demonstrate  reconstruction 
of  an  ideal  bandlimited  deterministic  tests  signal  and  a  bandlimited  random  process. 


Research  into  representing  a  function  by  its  sample  values,  and  development  of  the 
corresponding  interpolation  formulas,  enjoys  a  rich  history  [1-11].  Representations  may 
be  chosen,  e.g.  by  use  of  Fourier  series  coefficients,  that  do  not  restrict  the  frequency 
domain  support.  Although  Fourier  coefficients  are  attractive  for  many  reasons,  a  variety 
of  interpolation  functions  have  been  invoked  [1,3, 4, 7].  Other  convenient  representations, 
e.g.  use  of  the  function  sample  values  directly,  require  the  function  to  be  bandlimited. 
Bandlimited  functions  that  are  sampled  at  less  than  the  Nyquist  rate  [2]  exhibit  a 
distortion  termed  aliasing,  however,  many  situations  allow  the  aliasing  to  be  eliminated 
or  reduced  if  additional  information  is  available  at  the  sampling  instants  [8,9,1 1,12,15- 
18]. 

Staggered  Sampling  Theory 

Although  so  called  ‘Staggered  Sampling  Theory,”  [c.f.  20-35  as  a  small  reference 
sample]  is  outside  the  scope  of  this  study,  we  note  the  parallels  to  our  work.  Staggered 
sampling  is  widely  used  in  the  high-  speed  oscilloscope  industry  so  that  two  or  more 
parallel  streams  of  sampled  data,  offset  in  time,  may  be  used  to  lower  the  system¬ 
sampling  rate. 

Staggered  sampling  is  useful  for  several  applications  such  as  sensor  networks,  analog-to- 
digital  converter  (ADC),  and  sampling  high  bandwidth  analog  signals. 

Staggered  sampled  ADCs  are  an  alternative  for  parallelizing  the  sampling  task  over 
several  converters.  Such  a  system  may  provide  a  relatively  high  sampling  rate  with 
component  ADCs  operating  at  lower  rates.  The  issues  with  staggered  sampling  tend  to  be 
the  difficulty  of  maintaining  close  timing  accuracy  between  the  sample  streams  and 
quantization  errors  from  finite  arithmetic.  Conventional  staggered  sampling  theory  makes 
use  of  filter  bank  theory  with  required  up  and  down  sampling  converters  that  may  be  a 
disadvantage. 

In  this  research  we  provide  an  alternative  formulation  that  admits  sampling  directly  at  Vi 
the  Nyquist  rate  without  any  required  up  and  down  sampling.  Although  we  restrict  our 
attention  to  the  two-channel  case,  there  is  no  mathematical  reason  the  theory  cannot  be 
extended  to  a  greater  number  of  channels  (for  example  we  could  sample  at  Va  the  Nyquist 
rate  and  carry  four  closely  spaced  samples  each  system  sampling  instant). 

Impulse-doublet  Sampling 

A  generalised  version  of  the  Nyquist  sampling  theorem  [8]  admits  sampling  at  an  average 
rate  equal  to  twice  the  highest  frequency  component  of  the  sampled  signal.  For  example, 
we  may  envision  a  sampling  structure  wherein  a  pair  of  closely  spaced  impulses  perform 
the  signal  sampling,  each  impulse  pair  sampling  the  signal  at  one-half  the  minimum  rate 
required  for  ideal  impulse  sampling.  Let  f(t)  be  a  low  pass  signal  with  bandwidth  fc  =  W 
Hz  that  is  sampled  by  the  function 


p(t)  =  5(t  +  x/2)  +  5(t  -  x/2) 


at  a  rate  f s  =  f c  =  1/T  Hz,  0  <  x  «  T/2.  We  assume  that  f(t)  is  real  so  that  the  sampled 

spectrum  is  Hermitian,  and  may  be  restored  using  positive  frequencies  only.  The  aliased 

t  ....... 

spectrum,  F  (CO)  over  0  <  f  <  fc  (using  co  =  27tf  to  avoid  notational  difficulties)  is  given  by 

F*(co)  =  A0  F(oo)  +  AiF(co  -  cos) 


where 


Am  =  (2/T)  cos(mcocx/2) 

as  shown  in  Figure  1 . 

Reconstruction  Equations 

We  next  derive  a  frequency  domain  restoration  algorithm  and  its  corresponding  time 
domain  interpolation  formula. 

The  key  mathematical  manipulations  will  be  shown. 

F*(co)  =  Ao  F(co)  +  AiF(co  -  cos) ,  0  <  co  <  cocLet  coo  =  coc/2 

F*(coo)  =  Ao  F(coo)  +  AiF(coo  -  oos) 

Assume  f(t)  is  real  valued,  then 

F(-co)  =  Fa(co) 

F*(coo)  =  Ao  F(coo)  +  AiFa(cos  -  coo) 

But  cos  =  2  coo 

Therefore, 

F*(coo)  =  Ao  F(COO)  +  AiFA(coo) 

F*(coO)  =  Fr*(coO)  +  jFi*(coO) 

and 

F(coQ)  =  Fr(coO)  +  jFi  (coQ) 


It  follows  that 


Fr*(coO)  +  jFi*(®0)  =  Ao[Fr(coO)  +  jFi  (©0)]  +  Al[Fr(coO)  -  jFi  (coO)] 

Equating  the  real  and  imaginary  parts,  we  have 

Fr*(coO)  =  A()Fr(coO)  +  Ai  Fr  (coo) 

Fi*(coO)  =  AoFi(coO)  -  Al  Fi  (coO) 

Now,  consider 

co  =  coO  +/-  Aco 

in  the  aliased  overlap  region 

F*(coo  +/-  Aco)  =  Ao  F(coO  +/-  Aco)  +  AiF(coO  +/-  Aco  -  cos), 

0  <  co  <  (Oc 
Let 

coo  =  cOc/2  or  cos  =  2coO 

F*(coo  +/-  Aco)  =  Ao  F(coo+/-  Aco)  +  AiF(-coo  +/-  Aco) 
Assume  f(t)  is  real  valued,  then 

F(-co)  =  FA(C0) 

and 

F*(tt>0  +/-  Aco)  =  Ao  F(coo+/-  Aco)  +  AlFA(coo  -/+  Aco) 
Equating  Real  and  Imaginary  parts  results  in 

Fr*(coo  +/-  Aco)  =  Ao  Fr(coo+/-  Aco)  +  AiFrA(coo  -/+  Aco) 
Fi*(coo  +/-  Aco)  =  Ao  Fi(a>0+/-  Aco)  -  AiFiA(coo  -/+  Aco) 

Real  part: 

Fr*(coo  +  Aco)  =  Ao  Fr(coo+  Aco)  +  AlFrA(coo  -  Aco) 

Fr*(coo  -  Aco)  =  Ao  Fr(C0Q-  Aco)  +  AlFrA(coo  +  Aco) 


Imaginary  part: 


Fi*(coo  +  Aco)  =  Ao  Fi((00+  Aco)  -  AlFiA(a>0  -  Aco) 

Fi*(ooo  -  Aco)  =  Ao  Fi(a>0-  Aco)  -  AlFiA(coo  +  Aco) 

Solution  of  the  simultaneous  equations  leads  to 

Fr(a>0  +/-  Aco)  =  [Ao  Fr*(coo+/-  Aco)  -  AiFr*(coo  -/+  Aco)  ]/  [Ao^  -Al2] 

Fi(a>0  +/-  Aco)  =  [Ao  Fi*(coo+/-  Aco)  +  AiFi*(ooo  -/+  Aco)]  /  [Ao^  -Al2] 

Let  Aco  be  some  arbitrary  frequency  offset  from  co  0 

F*(co  o  +  Aco)  =  Ao  F(co  o  +  Aco)  +  AiF(co  o  -  Aco) 

F*(co  o  -  Aco)  =  Ao  F(co  o  -  Aco)  +  AiF(co  o  +  Aco) 

The  derivation,  by  straightforward  manipulation  of  the  aliased  spectral  equations  as  just 
derived  manually,  may  be  summarized  compactly  as 

F  =  A_1F* 

with 


A  = 


A 

4 


F(a>0  +  A  co) 
F(co0  -  A  co) 


F* 


F*(a)0  +  A  co) 
F  (O)0  -  A co) 


where  (Oo  =  ooc/2,  and  Aco  is  an  arbitrary  frequency  offset  from  (Oq. 


And  re-substituting 


(Os  =  2  coo 


F(co)  =  Ao/(Ao2  -  Ai2)  F*(co)  -  A\/  (Ao2  -  Ai2)  F*(co  -  cos) 

0  <  co  <  co  c 

Upon  application  of  a  brick  wall  ideal  filter,  G(co)  =  l,0<co<coC5 
to  the  positive  frequency  components  defined  by  the  frequency  domain  reconstruction 
formula,  the  time  domain  interpolation  equation  follows  directly  as 


F(co)  =  Ao/(Ao2  -  Ai2)  F*(co)  -  Ai /  (Ao2  -  A(2)  F*(co  -  cos) 

0  <  co  <  co  c 


f(t)  =  2Re { g(t) * ( A0/ ( A02  -  A(2)  f*(t)  -  A\!  (A02  -  Ai2)  f*(t)  e  j  “sty]} 


where  g(t)  =  F'l  (G(co)  } 
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where  f  (t)  is  the  sampled  version  of  f(t). 

Alternatively,  the  aliasing  could  be  removed  in  the  discrete  frequency  domain.  A 
frequency  domain  reconstruction  might  be  attractive  computationally,  however,  careful 
consideration  of  windowing  effects  would  be  required.  Another  advantage  may  be  that 
diagonal  loading  of  the  matrix  equations  could  lead  to  a  Maximum  Likelihood  solution  to 
mitigate  jitter  noise  effects.  Such  topics  are  outside  the  scope  of  this  investigation. 


We  may  quantify  the  effect  of  impulse  pair  spacing  on  the  numerical  conditioning  of  the 
restoration  algorithm.  The  eigenvalues  of  matrix  A  are  Ao  +/-  Ai,  Ao,  A]  >  0,  Ai  <  Ao; 
therefore  the  matrix  Condition  Number  of  A  is  given  by 


C.N. 


\  +  A 

A0  —  A1 


A 


o 


2 
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4  = 


—cos (tAVt)  . 
T 


We  observe  that: 

•  As  T  approaches  zero  the  C.N.  — >  as  we  expect,  i.e.,  matrix  A  is  singular,  and 

no  reconstruction  is  possible 

•  As  T  approaches  T/2  the  C.N.  — »  1,  i.e.,  Ai  =  0,  the  aliasing  is  zero,  and 
conventional  ideal  impulse  sampling  at  the  Nyquist  rate  obtains  in  the  limit 

•  For  0  <  x  <  T/2  numerical  stability  of  the  time  domain  interpolation  formula  or  the 
frequency  domain  alias  removal  algorithms  are  clearly  a  function  of  the  impulse- 
doublet  spacing. 

Simulation  Results 

The  following  simulations  illustrate  that  reconstruction  is  possible  with  the  aliased 
components  removed.  We  consider  an  ideal  bandlimited  test  signal,  the  “Sine  Function,” 
as  well  as  a  more  realistic  signal  modelling  a  bandlimited  random  process. 

Figure  2.  illustrates  the  results  of  sampling  a  .1  FIz  bandwidth  “Sine  Function,”  at  .1  FIz. 
We  reconstruct  100  points,  with  each  reconstructed  point  estimated  from  a  500-point 
summation.  The  impulse  doublet  has  width  .01  Second  for  this  example.  No  noise  has 
been  added. 

The  next  two  examples  illustrates  that  reconstruction  of  a  random  process  appears 
possible  (indeed  we  have  been  able  to  prove  convergence).  White  Gaussian  Noise  was 
filtered  to  .1  FIz  with  a  500  tap  FIR  filter.  The  sampling  rate  is  .1  Hz.  As  for  the  previous 
example,  the  impulse  doublet  width  is  set  arbitrarily  to  .01  second.  In  Figure  3.  we 
illustrate  the  reconstruction  of  this  random  process.  Figure  4.  shows  reconstruction  of  a 
noise-like  waveform,  16-QAM  with  added  noise.  Although  not  shown,  reconstruction  of 
16-QAM  without  added  noise  is  error  free.  In  fact  all  waveforms  tried,  even  a  pure 


sinusoid  (bandpass  rather  than  low  pass)  reconstructed  error  free  without  added  noise, 
somewhat  a  surprise  as  we  assumed  low  pass  spectra  during  the  derivation. 

The  next  simulations  demonstrate  the  effects  of  quantization  noise  and  timing  jitter 
between  the  staggered  samples  represented  by  the  doublet  impulse  pair.  We  also  illustrate 
the  effects  of  reconstruction  inaccuracies  due  to  quantization  on  an  ambiguity  function, 
widely  used  in  radar  signal  design  [191,  computed  from  a  Sine  function. 

Figures  5.-9.  illustrate  the  degrading  effects  on  computing  an  ambiguity  function  of  the 
Sine  function  for  8  bit,  16  bit,  and  32  bit  quantization.  The  Sine  function  is  chosen  as  an 
ideal  low  pass  waveform  that  possessing  infinite  time-domain  support,  thus  it’s  a  difficult 
signal  to  reconstruct.  The  results  are  not  surprising  although  might  have  desired  less 
reconstruction  error  for  16  bit  arithmetic.  The  32  bit  quantization  results  in  a  zero  error 
reconstruction. 

It  is  easier  to  view  the  degradation  due  to  jitter  and  quantization  direction  on  the 
reconstructed  Sine  function  rather  than  adding  the  extra  transformation  into  the  ambiguity 
plane.  Therefore,  the  remaining  simulations  focus  on  the  time  domain  reconstructed  Sine 
function  and  several  example  Power  Spectral  Densities  (PSDs)  to  show  the  impact  in 
frequency. 

For  the  results  in  Figures  10.  -  18.  we  fix  the  quantization  at  16  bits,  and  then  allow  the 
spacing  between  the  impulse  pair  to  vary.  Recall  that  we  showed  the  numerical  ill 
conditioning  directly  related  to  impulse  pair  spacing  with  reduced  spacing  resulting  in 
greater  ill  conditioning. 

A  spacing  of  .001  seconds  is  disastrous,  .005  seconds  much  improved,  and  degradation 
just  noticeable  at  a  spacing  of  .01  seconds.  At  a  spacing  of  .025  -  1.0  seconds  the 
reconstruction  error  is  zero.  Whether  this  degradation  is  acceptable  would  require 
analysis  and  simulation  based  upon  a  specific  radar  system. 

Fortunately,  from  Figures  19.  -  21.,  we  see  that  switching  to  32  bit  quantization  removes 
any  error  even  at  the  closely  spaced  case  of  .001  seconds. 

Less  fortunately,  although  unsurprising,  in  Figures  22.  -  26.  we  note  that  the  results  for  8 
bit  quantization  are  unacceptable. 

Figures  27.  and  28.  show  the  effect  of  jitter  with  a  doublet  spacing  of  .01  second,  and 
jitter  variances  of  .001  and  .005  with  predictable  results.  The  variance  level  was  chosen 
to  illustrate  an  acceptable  degradation  and  unacceptable  degradation  rather  than  any 
connection  to  a  technical  specification.  Again,  whether  these  jitter  levels  are  acceptable 
or  practical  requires  simulation  based  upon  real-world  hardware  specifications. 

Finally,  in  Figures  29.  -  31.  we  show  the  effects  of  jitter  on  the  PSD  of  the  reconstructed 
Sine  function.  Although  passband  distortion  is  apparent  the  major  component  of  the 
effects  of  jitter  appear  to  be  in  the  level  of  the  out  of  band  spectral  energy.  This  result 


bears  some  future  caution  and  study  into  the  effects  of  timing  jitter  on  the  reconstruction 
equation  accuracy. 


Conclusions  and  Further  Study 

We  have  summarised  time  domain  and  frequency  domain  restoration  algorithms  for  a 
signal  sampled  by  a  periodic  impulse-doublet  at  a  rate  equal  to  one-half  the  conventional 
Nyquist  rate.  The  numerical  stability  of  the  proposed  solution  is  a  function  of  the  C.N.  of 
matrix  A,  directly  related  to  the  impulse-doublet  spacing.  The  reconstruction  equations 
become  ill  posed  as  x  — >  0.  We  presented  simulation  results  that  demonstrate  the  accuracy 
of  the  time  domain  reconstruction  equations  with  and  without  noise  added  to  the  samples, 
time  domain  jitter  between  the  staggered  samples,  and  with  quantization  noise  added.  We 
note  that  there  does  not  appear  to  be  any  mathematical  reason  the  sampling  rate  cannot  be 
reduced,  and  the  additional  aliased  spectra  restored  in  an  extended  solution  to  the 
simplified  case  we  considered  here  of  first  order  aliasing. 

The  aliasing  could  also  be  removed  in  the  discrete  frequency  domain.  A  frequency 
domain  reconstruction  might  be  attractive  computationally,  however,  careful 
consideration  of  windowing  effects  would  be  required.  Another  advantage  may  be  that 
diagonal  loading  of  the  matrix  equations  could  lead  to  a  Maximum  Likelihood  solution  to 
mitigate  jitter  noise  effects.  Such  topics  are  outside  the  scope  of  this  investigation. 

Based  upon  the  simulation  results  applied  to  reconstructing  a  bandlimited  random 
process,  we  came  to  believe  that  the  restoration  equations  presented  here  are  indeed  valid 
for  bandlimited  random  processes.  Indeed,  we  have  demonstrated  mathematically  in  this 
research  that  reconstruction  of  an  aliased  bandlimited  random  process  is  possible  in  a 
“Limit  in  the  Mean”  sense. 

The  effect  of  A/D  quantisation  noise,  and  sample  timing  jitter,  should  be  quantified 
objectively  via  simulation  for  specific  hardware  and  systems,  before  the  potential  use  of 
these  reconstruction  algorithms  can  be  considered  for  an  operational  radar  system.  For 
off-line  waveform  design  purposes,  noise  sensitivities  may  be  less  of  an  issue. 
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Figure  1.  Aliased  Spectrum  Sampled  at  One-Half  the  Nyquist  Rate 
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Figure  2.  Reconstructed  “Sine  Function,”  500  Samples  Used,  x  =  .01  Sec,  No  Added 
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Band  Limited  Noise,  W=.1  Hz,  Fs=.1  Hz,  (Solid  —  Original,  Dots  —  Reconstructed) 


Figure  3. 

Figure  3.  Reconstructed  Bandlimited  Noise,  x  =  .01  Second 
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Figure  4.  Reconstructed  16-QAM  With  Additive  Noise,  x  =  2  Seconds,  Noise  Variance 
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Figure  5.  Ambiguity  Function  of  Sine  Function,  No  Noise 
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Figure  6.  Ambiguity  Function  of  Reconstructed  Sine  Function,  No  Noise 
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Figure  7.  Ambiguity  Function  of  Reconstructed  Sine  Function,  8  Bit  Quantization,  tau 

1.0 


Ambiguity  Function  of  Reconstructed  Signal 


0.25 

0.2 

0.15 

0.1 

0.05 

0 

250 


250 


Figure  8.  Ambiguity  Function  of  Reconstructed  Sine  Function,  16  Bit  Quantization 
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Figure  9.  Ambiguity  Function  of  Sine  Function,  32  Bit  Quantization 
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Figure  10.  Reconstructed  Sine  Function,  16  Bit  Quantization 
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Figure  11.  Reconstructed  Sine  Function,  16  Bit  Quantization!,  Tau  =  .001  Sec 
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Figure  12.  Reconstructed  Sine  Function,  16  Bit  Quantization,  Tau  =  .005  Sec 
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Figure  13.  Reconstructed  Sine  Function,  16  Bit  Quantization,  Tau  =  .005  Sec 
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Figure  14.  Reconstructed  Sine  Function,  16  Bit  Quantization,  Tau  =  .01  Sec 
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Figure  15.  Reconstructed  Sine  Function,  16  Bit  Quantization,  Tau  =  .025  Sec 
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Figure  16.  Reconstructed  Sine  Function,  16  Bit  Quantization,  Tau  =  .025  Sec 
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Figure  17.  Reconstructed  Sine  Function,  16  Bit  Quantization,  Tau  =  .05  Sec 
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Figure  18.  Reconstructed  Sine  Function,  16  Bit  Quantization,  Tau  =  1.0  Sec 
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Figure  19.  Reconstructed  Sine  Function,  32  Bits  Quantization,  Tau  =  .001  Sec 
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Figure  20.  Reconstructed  Sine  Function,  32  Bit  Quantization,  Tau  =  .005  Sec 
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Figure  21.  Reconstructed  Sine  Function,  32  Bit  Quantization,  Tau  =  .01  Sec 
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Figure  22.  Reconstructed  Sine  Function,  8  Bit  Quantization,  Tau  =  .05  Sec 
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Figure  23.  Reconstructed  Sine  Function,  8  Bit  Quantization,  Tau  =  .01  Sec 
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Figure  24.  Reconstructed  Sine  Function,  8  Bit  Quantization,  Tau  =  .025  Sec 
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Figure  25.  Reconstructed  Sine  Function,  8  Bit  Quantization,  Tau  =  .05  Sec 
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Figure  26.  Reconstructed  Sine  Function,  8  Bit  Quantization,  Tau  =1.0  Sec 
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Figure  27.  Reconstructed  Sine  Function,  Jitter  Variance  =  .001,  Tau  =  .01  Sec 
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Figure  28.  Reconstructed  Sine  Function,  Jitter  Variance  =  .005,  Tau  =  .01  Sec 
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Figure  29.  PSD  of  Sine,  Jitter  Variance  =  .001,  Tau  =  .01  Sec 
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Figure  30.  PSD  of  Sine,  Jitter  Variance  =  .005,  Tau  =  .01  Sec 


Figure  31.  PSD  of  Sine,  No  Jitter,  Tau  =  .01  Sec 
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1Abstract  -  A  primary  goal  in  radar 
waveform  design  is  to  achieve  a 
prescribed  discrete  time  radar 
ambiguity  function.  Radar  waveform 
design  to  meet  a  given  ambiguity 
function  specification  is  an  important 
issue  under  the  waveform  diversity 
technology.  Much  of  the  development 
in  this  area  is  for  continuous  time 
signals.  For  the  discrete  time  case,  the 
sampling  rate  becomes  an  important 
design  issue.  On  the  one  hand,  too  low 
of  a  sampling  rate  results  in  spectral 
aliasing,  on  the  other  hand,  a 
sampling  rate  chosen  higher  than 
necessary  increases  the  computational 
burden.  We  show  in  this  paper  that 
aliased  spectra,  arising  from  sampling 
below  the  Nyquist  rate,  may  be 
completely  eliminated. 


1  The  USAF  AFOSR  AOARD,  under  Contract  # 
044046,  “Radar  Waveform  Design  of 
Undersampled  Signals”  is  supporting  this 
research. 


I.  INTRODUCTION 
For  the  purpose  of  this  research,  we 
assume  for  a  start  that  the  desired 
waveform  is  available.  Since  radar 
processing  is  done  digitally,  it  is 
important  to  examine  the  effects  of 
sampling  the  ambiguity  function  in  the 
Delay-Doppler  plane  and  study  the 
resulting  resolution  trade-offs, 
reconstruction  issues,  and  aliasing 
problems.  The  preliminary  techniques 
developed  in  this  research,  and  applied 
to  the  ambiguity  function,  should  be 
quite  useful  in  this  context. 

The  approach  for  this  research  will  be  to 
apply  the  proposed  aliased  signal 
reconstruction  algorithms  to  radar 
waveform  design  against  ambiguity 
plane  constraints.  We  hypothesise,  based 
on  simulation  results,  that  the 
reconstruction  algorithm  is  valid  for 
bandlimited  random  processes. 

Another  component  of  the  research  is  to 
prove  analytically  that  exact 
reconstruction  of  bandlimited  random 


processes  is  possible  in  a  “Limit  in  the 
Mean”  sense.  This  is  an  important  result 
since  modem  radar  waveform  design 
involves  the  use  of  pseudorandom 
coding  for  pulse  compression  and 
LPI/LPD  waveforms. 

In  addition  to  analytically  proving 
convergence  of  the  reconstruction 
algorithm,  we  plan  to  present  design 
trade-off  plots  using  Matlab-based 
simulation  tools. 

The  effect  of  A/D  quantisation  noise, 
impulse  doublet  spacing,  and  sample 
timing  jitter,  need  to  be  quantified  via 
simulation,  before  the  potential  use  of 
these  reconstruction  algorithms  can  be 
considered  for  an  operational  radar 
system.  Matlab-based  simulation  tools 
will  be  utilised  to  develop  reconstruction 
sensitivities  to  the  various  real-world 
degradations  listed  above.  For  off-line 
waveform  design  purposes,  noise 
sensitivities  may  be  less  of  an  issue. 

II.  APPLICATIONS  TO  RADAR 
WAVEFORM  DESIGN 
A  primary  goal  in  radar  waveform 
design  is  to  achieve  a  prescribed  discrete 
time  radar  ambiguity  function.  Radar 
waveform  design  to  meet  a  given 
ambiguity  function  specification  is  an 
important  issue  under  the  waveform 
diversity  technology.  Much  of  the 
development  in  this  area  is  for 
continuous  time  signals.  For  the  discrete 
time  case,  the  sampling  rate  becomes  an 
important  design  issue.  On  the  one  hand, 
too  low  of  a  sampling  rate  results  in 
spectral  aliasing,  on  the  other  hand,  a 
sampling  rate  chosen  higher  than 
necessary  increases  the  computational 
burden.  We  show  in  this  research  project 
that  aliased  spectra,  arising  from 


sampling  below  the  Nyquist  rate,  may  be 
completely  eliminated. 

The  radar  ambiguity  function  based 
waveform  design  is  complicated  by  the 
fact  that  for  a  given  waveform,  the 
ambiguity  function  can  be  readily 
calculated.  However,  given  an  ambiguity 
function  specification,  it  is  possible  to 
have  more  than  one  waveform  that  meets 
the  specification.  For  the  purpose  of  this 
research,  we  assume  for  a  start  that  the 
desired  waveform  is  available.  Since 
radar  processing  is  done  digitally,  it  is 
important  to  examine  the  effects  of 
sampling  the  ambiguity  function  in  the 
Delay-Doppler  plane  and  study  the 
resulting  resolution  trade-offs, 
reconstruction  issues,  and  aliasing 
problems.  The  preliminary  techniques 
developed  in  this  research,  and  applied 
to  the  ambiguity  function  in  this 
proposed  research,  should  be  quite 
useful  in  this  context. 

The  goal  of  this  research  is  to  apply  the 
aliased  signal  reconstruction  algorithm 
to  radar  waveform  design  against 
ambiguity  plane  constraints  [19].  The 
radar  receiver  discrete  time  matched 
filter  computational  complexity  may  be 
potential  reduced  by  implementing  the 
signal  restoration  algorithm  summarised 
in  this  proposal.  We  hypothesise,  based 
on  simulation  results  presented  here,  that 
the  reconstruction  algorithm  is  valid  for 
bandlimited  random  processes. 

III.  SAMPLING  THEORIES 
Research  into  representing  a  function  by 
its  sample  values,  and  development  of 
the  corresponding  interpolation 
formulas,  enjoys  a  rich  history  [1-11]. 
Representations  may  be  chosen,  e.g.  by 
use  of  Fourier  series  coefficients,  that  do 
not  restrict  the  frequency  domain 


support.  Although  Fourier  coefficients 
are  attractive  for  many  reasons,  a  variety 
of  interpolation  functions  have  been 
invoked  [1,3, 4, 7].  Other  convenient 
representations,  e.g.  use  of  the  function 
sample  values  directly,  require  the 
function  to  be  bandlimited.  Bandlimited 
functions  that  are  sampled  at  less  than 
the  Nyquist  rate  [2]  exhibit  a  distortion 
termed  aliasing,  however,  many 
situations  allow  the  aliasing  to  be 
eliminated  or  reduced  if  additional 
information  is  available  at  the  sampling 
instants  [8,9,11,12,15-18]. 

IV.  IMPULSE-DOUBLET  SAMPLING 
A  generalised  version  of  the  Nyquist 
sampling  theorem  [8]  admits  sampling  at 
an  average  rate  equal  to  twice  the 
highest  frequency  component  of  the 
sampled  signal.  For  example,  we  may 
envision  a  sampling  structure  wherein  a 
pair  of  closely  spaced  impulses  perform 
the  signal  sampling,  each  impulse  pair 
sampling  the  signal  at  one-half  the 
minimum  rate  required  for  ideal  impulse 
sampling.  Let  f(t)  be  a  low  pass  signal 
with  bandwidth  fc  =  W  Hz  that  is 
sampled  by  the  function 

p(t)  =  8(t  +  X/2)  +  8(t  -  X/2) 

at  a  rate  fs  =  fc  =  1/T  Hz,  0  <  x  «  T/2. 
We  assume  that  f(t)  is  real  so  that  the 
sampled  spectrum  is  Hermitian,  and  may 
be  restored  using  positive  frequencies 
only.  The  aliased  spectrum,  F  (CO)  over  0 
<  f  <  fc  (using  (0  =  2nf  to  avoid 
notational  difficulties)  is  given  by 

F*((0)  =  A0  F(co)  +  AiF(co  -  C0S) 
where 


As  shown  in  Figure  1 . 

V.  RECONSTRUCTION  EQUATIONS 
We  next  derive  a  frequency  domain 
restoration  algorithm  and  its 
corresponding  time  domain  interpolation 
formula.  Straightforward  manipulation 
of  the  aliased  spectral  equations  results 
in  the  solution 


F  =  A_1F* 


with 


A  = 


A0  A{ 

.A  A 


F(co0  +  Ary) 
F(cq0  -  A co) 


F*  = 


F\co0  +  Ary) 
F\o) o  -  Ary) 


where  Oio  =  (0c/2,  and  Aco  is  an  arbitrary 
frequency  offset  from  COo. 

The  interpolation  equation  follows 
directly  as 


/(f)  =  2  Re- 


x{nT  +  r  /  2)g(t  —  nT  +  z  12) 

+  x(nT -z/2)g(t-nT  -z/2) 


with 


g(0  = 


W  sin(^Wt)  m, 
xWt 


Am  =  (2/T)  cos(mcocx/2) 


x(t)  = 


A0-Aiei2mt 
A2  -  A2 


fit) 


where  f  (t)  is  the  sampled  version  of  f(t). 


Finally,  we  quantify  the  effect  of 
impulse  pair  spacing  on  the  numerical 
conditioning  of  the  restoration 
algorithm. 


VI.  SIMULATION  RESULTS 
The  following  simulations  illustrate  that 
reconstruction  is  possible  with  the 
aliased  components  removed.  We 
consider  an  ideal  test  signal,  the  “Sine 
Function,”  as  well  as  a  realistic  signal, 
Bandlimited  Noise,  and  its  ambiguity 
function.  We  also  show  the  ambiguity 
function  of  the  reconstructed  sine 
function. 


The  eigenvalues  of  matrix  A  are  Ao  +/- 
Ai,  Ao,  Ai  >  0;  therefore  the  matrix 
Condition  Number  of  A  is  given  by 


C.N. 


A)  +  A 
A) —  A 


A 


o 


2_ 

T 


A  = 


—  cos (7tWt)  . 
T 


Figure  2.  shows  the  results  of  sampling  a 
.1  Hz  bandwidth  “Sine  Function,”  at  .1 
Hz.  We  reconstruct  100  points,  with 
each  reconstructed  point  estimated  from 
a  500  point  summation.  The  impulse 
doublet  has  width  .1  Second  for  this 
example.  No  noise  has  been  added. 

The  next  example  illustrates  that 
reconstruction  of  a  random  process  is 
possible.  White  Gaussian  Noise  was 
filtered  to  .1  Hz  with  a  500  tap  FIR 
filter.  The  sampling  rate  is  .1  Hz.  In 
Figure  3.  we  illustrate  the  reconstruction 
of  this  random  process. 


We  observe  that: 

•  As  T  approaches  zero  the  C.N.  — > 
°°  as  we  expect,  i.e.,  no 
reconstruction  is  possible 


Figures  4.,  5.,  and  6.,  show  the 
equivalence  of  the  ambiguity  functions 
computed  from  the  original  signals  or 
reconstructed  signals  from  aliased 
copies. 


•  As  T  approaches  T/2  the  C.N.  — > 
1,  i.e.,  Ai  =  0,  the  aliasing  is 
zero,  and  conventional  ideal 
impulse  sampling  at  the  Nyquist 
rate  obtains  in  the  limit 

•  For  0  <  T  <  T/2  numerical 
stability  of  the  time  domain 
interpolation  formula  or  the 
frequency  domain  alias  removal 
algorithms  are  clearly  a  function 
of  the  impulse-doublet  spacing. 


VII.  CONCFUSIONS 
We  have  summarised  time  domain  and 
frequency  domain  restoration  algorithms 
for  a  signal  sampled  by  a  periodic 
impulse-doublet  at  a  rate  equal  to  one- 
half  the  conventional  Nyquist  rate.  The 
numerical  stability  of  the  proposed 
solution  is  a  function  of  the  C.N.  of 
matrix  A,  directly  related  to  the  impulse- 
doublet  spacing.  The  problem  becomes 
ill  posed  as  x  —>  0.  We  presented 
simulation  results  that  demonstrate  the 
accuracy  of  the  time  domain 
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reconstruction  equations  with  and 
without  noise.  Finally,  we  note  that  there 
does  not  appear  to  be  any  mathematical 
reason  the  sampling  rate  cannot  be 
reduced,  and  the  additional  aliased 
spectra  restored  in  an  extended  solution 
to  the  simplified  case  we  considered  here 
of  first  order  aliasing. 

Based  upon  the  simulation  results 
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The  effect  of  A/D  quantisation  noise, 
and  sample  timing  jitter,  need  to  be 
quantified  via  simulation,  before  the 
potential  use  of  these  reconstruction 
algorithms  can  be  considered  for  an 
operational  radar  system.  For  off-line 
waveform  design  purposes,  noise 
sensitivities  may  be  less  of  an  issue. 

We  propose  to  extend  these  results  to 
reducing  or  eliminating  the  aliasing  in  an 
ambiguity  function  computed  from  an 
undersampled  radar  waveform.  This 
proposed  research  would  lead  to  a  better 
understanding  of  the  effects  of  sampling, 
and  aliasing,  on  the  ambiguity  function 
computed  from  a  digitised  radar 
waveform. 
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Fig.  1.  Aliased  Spectrum  Sampled  at  One-Half  the  Nyquist  Rate 


f(t),  fhat(t) 


Sine  Eg,  W=.1  Hz,  Fs=.1  Hz,  (Solid  —  Original,  Dots  —  Reconstructed) 


Fig.  2.  Reconstructed  “Sine  Function,”  500  Samples  Used,  x  =  .01  Sec,  No  Added  Noise 


Band  Limited  Noise,  W=.1  Hz,  Fs=.1  Hz,  (Solid  —  Original,  Dots  —  Reconstructed) 


Fig.  3.  Reconstructed  Bandlimited  Noise,  x  =  .01  Second 
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Fig.  4.  Ambiguity  Function  of  Original  Signal  (.1  Hz  Bandlimited  Noise) 
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Fig.  5.  Ambiguity  Function  of  Reconstructed  Signal  (.1  Hz  Bandlimited  Noise) 
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Fig.  6.  Ambiguity  Function  of  the  “Sine”  Function 


Appendix  C:  Matlab  Code 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  .25  SQRTCOS  filtered  16-QAM 
%  BW  =  .2  Hz,  Sampled  at  .2  Hz 

elf 

clear 


Fs=l 

N=4097 

W=input( 'Enter  W:  ') 

T=fix(l/W) 

%tau=input('Enter  tau:  ') 
tau=2 

noise_var=input('Enter  noise  variance:  ') 
white_noise=sqrt(noise_var)  *randn(  1  ,N) ; 

AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
t=-50:50; 


load  -MAT  snr_inf.bak 
qam_sig=zeros(  1  ,N) ; 

qam_sig(2049-1250:2049+1250)=real(MatrixData(100:2600)); 

w=hamming(l,2501); 

%qam_sig(2049-1250:2049+1250)=w.*qam_sig(2049- 1250:2049+1250); 

f=zeros(l,N); 

xterml=zeros(l,N); 

xterm2=zeros(l,N); 

gterml=zeros(l,N); 

gterm2=zeros(l,N); 

f  =  qam_sig  +  white_noise; 

%  Construct  interpolation  functions  sampled  at  1  Hz 
for  n=l:Fs:2048 

xterml(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

xterm2(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

gterml(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

gterm2(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 


end 


for  n=2050:Fs:4097 

xterml(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

xterm2(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

gterml(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

gterm2(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

end 

xterm 1  ( 2049)=(  AO- A 1 )  *f(2049)/(A0A2- A 1 A2) ; 
xterm2(2049)=(  AO- A 1 )  *f(2049)/(A0A2- A 1 A2) ; 
gterml(2049)=W; 
gterm2 ( 2049)= W; 

%Interpolate  at  .2  Hz  sampling 
f_hat_save=zeros(l,length(t)); 
for  t_index=l :  length  (t) 
f_hat=zeros(l  ,length(t)); 

%for  m  =  T+length(t)+l:T:N-T-length(t) 
count=0; 

for  m  =  2049- 1250:T:2049+ 1250 

count=count+l; 

m_save(count)=m; 

f_hat(t_index)  =  f_hat(t_index)  +  xterml(m+tau/2)*gterml(-t(t_index)+m+tau/2)  + 
xterm2  ( m-tau/2)  *  gterm2  (-t(t_index)+m-tau/2) ; 

end 

f_hat_save(t_index)=2*real(f_hat(t_index)); 


end 


f_test=f(2049-fix(length(t)/2):2049+fix(length(t)/2)) 

f_hat_save 

snr=10*logl0(sum(qam_sig.A2)/noise_var) 

mse=sum((f_test-f_hat_save).A2) 

y=2*W*real(conv(xterml,gterml)  +  conv(xterm2,gterm2)); 
y(4097-fix(length(t)/2):4097+fix(length(t)/2)) 

plot(t,f_test) 
hold  on 

plot(t,f_hat_save,7) 

grid 

xlabel('Time  (Seconds)') 
ylabel('f(t),  fhat(t)') 

title('16-QAM  Eg,  W=.2  Hz,  Fs=.2  Hz,  (Solid  —  Original,  Dots  —  Reconstructed)') 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  Sine  function 
%  BW  =  .1  Hz,  Sampled  at  .1  Hz 


Fs=.l 
N=500 
Fc=.l 
W=.  1 


T=l/W 

tau=input('Enter  tau:  ') 
noise_var=input('Enter  noise  variance:  ') 
AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
t=-50:50 


for  t_index=l :  length  (t) 
t_index 

f=zeros(2,N); 
xterml=zeros(l,N); 
xterm2=zeros(l,N); 
gterml  =zeros(  1  ,N) ; 
gterm2=zeros(l,N); 


%  Construct  test  signal  and  interpolation  functions  sampled  at  1  Hz 

for  n=l:N 

arg=(n-N/2)*T; 

f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2))  + 
sqrt(noise_var)*randn(l,l); 

f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2))  + 
sqrt(noise_var)*randn(l,l); 

xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 
xterm2(n)=(A0-Al*exp(j*2*pi*W*(arg-tau/2)))*f(2,n)/(A0A2-AlA2); 
gterml  (n)=W*sin(pi*W*(t(t_index)-arg-tau/2))*exp(j*pi*W*(t(t_index)-arg- 
tau/2))/(pi*W*(t(t_index)-arg-tau/2)); 

gterm2(n)=W*sin(pi*W*(t(t_index)-arg+tau/2))*exp(j*pi*W*(t(t_index)- 

arg+tau/2))/(pi*W*(t(t_index)-arg+tau/2)); 


end 


%for  n=(N/2)+l:N 


%  arg=(n-249)*T; 

%  f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2)); 

%  f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2)); 

%  xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 

%  xterm2(n)=(  AO- A 1  *exp(j  *2  *pi*  W*  (arg-tau/2)))  *f(2,n)/(A0A2- A 1 A2) ; 

%  gterml(n)=W*sin(pi*W*(-arg+tau/2))*exp(j*pi*W*(-arg+tau/2))/(pi*W*(- 
arg+tau/2)); 

%  gterm2(n)=W*sin(pi*W*(-arg-tau/2))*exp(j*pi*W*(-arg-tau/2))/(pi*W*(-arg-tau/2)); 
%end 


%f(N/2)  =  2*W; 

%xterml  (N/2)=(  AO-  A 1 )  *f  (N/2)/(A0A2- A 1 A2) ; 
%xterm2(N/2)=(A0-Al)*f(N/2)/(A0A2-AlA2); 
%gterml(N/2)=W; 

%gterm2(N/2)=W; 


%  Interpolate  at .  1  Hz  sampling 

f_hat_0=0; 

for  m=l:N 

f_hat_0  =  f_hat_0  +  xterml(m)*gterml(m)  +  xterm2(m)*gterm2(m); 
end 

f_hat_0_save(t_index)=2*real(f_hat_0); 

%2*real(f_hat_0) 

%f(l,N/2) 

%f(2,N/2) 

end 

f_test=zeros(  1  ,length(t)) ; 


for  n=l:50 


f_test(n)  =  2*W*sin(2*pi*W*t(n))/(2*pi*W*t(n)); 
f_test(n+51)  =  2*W*sin(2*pi*W*t(n+51))/(2*pi*W*t(n+51)); 
end 

f_test(51)=.2; 

plot(t,f_test) 
hold  on 

plot(t,f_hat_0_save,'.') 

grid 

xlabel('Time  (Seconds)') 
ylabel('f(t),  fhat(t)') 

title('Sinc  Eg,  W=.l  Hz,  Fs=.l  Hz,  (Solid  —  Original,  Dots  —  Reconstructed)') 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  Sine  function 
%  BW  =  .1  Hz,  Sampled  at  .1  Hz 


Fs=.l 
N=500 
Fc=.l 
W=.  1 


T=l/W 

tau=input('Enter  tau:  ') 
noise_var=input('Enter  noise  variance:  ') 
AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
t=-50:50 


for  t_index=l :  length  (t) 
t_index 

f=zeros(2,N); 
xterml=zeros(l,N); 
xterm2=zeros(l,N); 
gterml  =zeros(  1  ,N) ; 
gterm2=zeros(l,N); 


%  Construct  test  signal  and  interpolation  functions  sampled  at  1  Hz 

for  n=l:N 

arg=(n-N/2)*T; 

f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2))  + 
sqrt(noise_var)*randn(l,l); 

f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2))  + 
sqrt(noise_var)*randn(l,l); 

xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 
xterm2(n)=(A0-Al*exp(j*2*pi*W*(arg-tau/2)))*f(2,n)/(A0A2-AlA2); 
gterml  (n)=W*sin(pi*W*(t(t_index)-arg-tau/2))*exp(j*pi*W*(t(t_index)-arg- 
tau/2))/(pi*W*(t(t_index)-arg-tau/2)); 

gterm2(n)=W*sin(pi*W*(t(t_index)-arg+tau/2))*exp(j*pi*W*(t(t_index)- 

arg+tau/2))/(pi*W*(t(t_index)-arg+tau/2)); 


end 


%for  n=(N/2)+l:N 


%  arg=(n-249)*T; 

%  f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2)); 

%  f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2)); 

%  xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 

%  xterm2(n)=(  AO- A 1  *exp(j  *2  *pi*  W*  (arg-tau/2)))  *f(2,n)/(A0A2- A 1 A2) ; 

%  gterml(n)=W*sin(pi*W*(-arg+tau/2))*exp(j*pi*W*(-arg+tau/2))/(pi*W*(- 
arg+tau/2)); 

%  gterm2(n)=W*sin(pi*W*(-arg-tau/2))*exp(j*pi*W*(-arg-tau/2))/(pi*W*(-arg-tau/2)); 
%end 


%f(N/2)  =  2*W; 

%xterml  (N/2)=(  AO-  A 1 )  *f  (N/2)/(A0A2- A 1 A2) ; 
%xterm2(N/2)=(A0-Al)*f(N/2)/(A0A2-AlA2); 
%gterml(N/2)=W; 

%gterm2(N/2)=W; 


%  Interpolate  at .  1  Hz  sampling 

f_hat_0=0; 

for  m=l:N 

f_hat_0  =  f_hat_0  +  xterml(m)*gterml(m)  +  xterm2(m)*gterm2(m); 
end 

f_hat_0_save(t_index)=2*real(f_hat_0); 

%2*real(f_hat_0) 

%f(l,N/2) 

%f(2,N/2) 

end 

f_test=zeros(  1  ,length(t)) ; 


for  n=l:50 


f_test(n)  =  2*W*sin(2*pi*W*t(n))/(2*pi*W*t(n)); 
f_test(n+51)  =  2*W*sin(2*pi*W*t(n+51))/(2*pi*W*t(n+51)); 
end 

f_test(51)=.2; 

plot(t,f_test) 
hold  on 

plot(t,f_hat_0_save,'.') 

grid 

xlabel('Time  (Seconds)') 
ylabel('f(t),  fhat(t)') 

title('Sinc  Eg,  W=.l  Hz,  Fs=.l  Hz,  (Solid  —  Original,  Dots  —  Reconstructed)') 
figure(2) 

amb_test=ambiguity(f_test,'true'); 

waterfall(amb_test) 

xlabel( 'T ime( +/- 100  Sec)') 

ylabel('Freq  (+/-  .25  Hz)') 

title('Ambiguity  Function  of  Original  Signal') 

figure(3) 

amb_test_hat=ambiguity(f_hat_0_save,'true'); 
waterfall(amb_test_hat) 
xlabel( 'T ime( +/- 100  Sec)') 
ylabel('Freq  (+/-  .25  Hz)') 

title('Ambiguity  Function  of  Reconstructed  Signal') 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  Sine  function  with  JITTER  on  sampling 
%  BW  =  .1  Hz,  Sampled  at  .1  Hz 


Fs=.l 
N=500 
Fc=.  1 
W=.l 
T=l/W 

tau=input('Enter  tau:  ') 
noise_var=input('Enter  noise  variance:  ') 
%b=input('Enter  number  of  bits:  ') 
%delta=2A-b; 

AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
t=-50:50 


for  t_index=l :  length  (t) 
t_index 

f=zeros(2,N); 
xterml=zeros(l,N); 
xterm2=zeros(l,N); 
gterml  =zeros(  1  ,N) ; 
gterm2=zeros(l,N); 


%  Construct  test  signal  and  interpolation  functions  sampled  at  1  Hz 

for  n=l:N 

arg=(n-N/2)*T; 

jit  1  =noise_var*randn(  1,1); 
jit2=noise_var*randn(  1,1); 

f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2  +  jitl))/(2*pi*W*(arg+tau/2  +  jitl)); 
f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2  +  jit2))/(2*pi*W*(arg-tau/2  +  jit2)); 

xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 
xterm2(n)=(A0-Al*exp(j*2*pi*W*(arg-tau/2)))*f(2,n)/(A0A2-AlA2); 
gterml  (n)=W*sin(pi*W*(t(t_index)-arg-tau/2))*exp(j*pi*W*(t(t_index)-arg- 
tau/2))/(pi*W*(t(t_index)-arg-tau/2)); 

gterm2(n)=W*sin(pi*W*(t(t_index)-arg+tau/2))*exp(j*pi*W*(t(t_index)- 

arg+tau/2))/(pi*W*(t(t_index)-arg+tau/2)); 


end 


%for  n=(N/2)+l:N 
%  arg=(n-249)*T; 

%  f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2)); 

%  f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2)); 

%  xterml  (n)=(  AO- A 1  *exp(j  *2  *pi*  W*  (arg+tau/2)))  *f  ( 1  ,n)/(A0A2- A 1 A2) ; 

%  xterm2(n)=(  AO- A 1  *exp(j  *2  *pi*  W*  (arg-tau/2)))  *f(2,n)/(A0A2- A 1 A2) ; 

%  gterml(n)=W*sin(pi*W*(-arg+tau/2))*exp(j*pi*W*(-arg+tau/2))/(pi*W*(- 
arg+tau/2)); 

%  gterm2(n)=W*sin(pi*W*(-arg-tau/2))*exp(j*pi*W*(-arg-tau/2))/(pi*W*(-arg-tau/2)); 
%end 


%f(N/2)  =  2*W; 

%xterml  (N/2)=(A0- A 1 )  *f  (N/2)/(A0A2- A 1 A2) ; 
%xterm2(N/2)=(A0-Al)*f(N/2)/(A0A2-AlA2); 
%gterml(N/2)=W; 

%gterm2(N/2)=W; 


%  Interpolate  at .  1  Hz  sampling 

f_hat_0=0; 

for  m=l:N 

f_hat_0  =  f_hat_0  +  xterml(m)*gterml(m)  +  xterm2(m)*gterm2(m); 
end 

fjiat_0_save(t_index)=2*real(f_hat_0); 

%2*real(f_hat_0) 

%f(l,N/2) 

%f(2,N/2) 


end 


f_test=zeros(  1  ,length(t)) ; 
for  n=l:50 

f_test(n)  =  2*W*sin(2*pi*W*t(n))/(2*pi*W*t(n)); 
f_test(n+51)  =  2*W*sin(2*pi*W*t(n+51))/(2*pi*W*t(n+51)); 
end 

f_test(51)=.2; 

figure(l) 
plot(t,f_test) 
hold  on 

plot(t,f_hat_0_save(l:101),'.') 

grid 

xlabel('Time  (Seconds)') 
ylabel('f(t),  fhat(t)') 

title('Sinc  Eg,  W=.l  Hz,  Fs=.l  Hz,  (Solid  —  Original,  Dots  —  Reconstructed)') 
hold  off 

figure(2) 

amb_test=ambiguity(f_test,'true'); 

waterfall(amb_test) 

xlabel( 'T ime( +/- 100  Sec)') 

ylabel('Freq  (+/-  .25  Hz)') 

title('Ambiguity  Function  of  Original  Signal') 

figure(3) 

amb_test_hat=ambiguity(f_hat_0_save,'trae'); 
waterfall(amb_test_hat) 
xlabel( 'T ime( +/- 100  Sec)') 
ylabel('Freq  (+/-  .25  Hz)') 

title('Ambiguity  Function  of  Reconstructed  Signal') 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  Sine  function 
%  BW  =  .1  Hz,  Sampled  at  .1  Hz 


Fs=.l 

N=500 

Fc=.l 

W=.075 

T=l/W 


tau=input('Enter  tau:  ') 
noise_var=input('Enter  noise  variance:  ') 
AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
t=-50:50 


for  t_index=l :  length  (t) 
t_index 

f=zeros(2,N); 
xterml=zeros(l,N); 
xterm2=zeros(l,N); 
gterml  =zeros(  1  ,N) ; 
gterm2=zeros(l,N); 


%  Construct  test  signal  and  interpolation  functions  sampled  at  1  Hz 

for  n=l:N 

arg=(n-N/2)*T; 

f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))  +  sqrt(noise_var)*randn(l,l); 
f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))  +  sqrt(noise_var)*randn(l,l); 

xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 
xterm2(n)=(A0-Al*exp(j*2*pi*W*(arg-tau/2)))*f(2,n)/(A0A2-AlA2); 
gterml  (n)=W*sin(pi*W*(t(t_index)-arg-tau/2))*exp(j*pi*W*(t(t_index)-arg- 
tau/2))/(pi*W*(t(t_index)-arg-tau/2)); 

gterm2(n)=W*sin(pi*W*(t(t_index)-arg+tau/2))*exp(j*pi*W*(t(t_index)- 

arg+tau/2))/(pi*W*(t(t_index)-arg+tau/2)); 

end 


%for  n=(N/2)+l:N 


%  arg=(n-249)*T; 

%  f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2)); 

%  f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2)); 

%  xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 

%  xterm2(n)=(  AO- A 1  *exp(j  *2  *pi*  W*  (arg-tau/2)))  *f(2,n)/(A0A2- A 1 A2) ; 

%  gterml(n)=W*sin(pi*W*(-arg+tau/2))*exp(j*pi*W*(-arg+tau/2))/(pi*W*(- 
arg+tau/2)); 

%  gterm2(n)=W*sin(pi*W*(-arg-tau/2))*exp(j*pi*W*(-arg-tau/2))/(pi*W*(-arg-tau/2)); 
%end 


%f(N/2)  =  2*W; 

%xtennl  (N/2)=(  AO- A 1 )  *f  (N/2)/(A0A2- A 1 A2) ; 
%xterm2(N/2)=(A0-Al)*f(N/2)/(A0A2-AlA2); 
%gterml(N/2)=W; 

%gterm2(N/2)=W; 


%  Interpolate  at .  1  Hz  sampling 

f_hat_0=0; 

for  m=l:N 

f_hat_0  =  f_hat_0  +  xterml(m)*gterml(m)  +  xterm2(m)*gterm2(m); 
end 

f_hat_0_save(t_index)=2*real(f_hat_0); 

%2*real(f_hat_0) 

%f(l,N/2) 

%f(2,N/2) 

end 

f_test=zeros(  1  ,length(t)) ; 
for  n=l:101 

f_test(n)  =  2*W*sin(2*pi*W*t(n)); 
end 


plot(t,f_test) 
hold  on 

plot(t,f_hat_0_save,'.') 

grid 

xlabel('Time  (Seconds)') 
ylabel('f(t),  fhat(t)') 

title('Sinewave  Eg,  W=.075  Hz,  Fs=.075  Hz,  (Solid  —  Original,  Dots  —  Reconstructed)') 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  Sine  function  —  Noise  Simulations 
%  BW  =  .1  Hz,  Sampled  at  .1  Hz 


Fs=.l 

N=500 

Fc=.  1 

W=.l 

T=l/W 

tau=.l 

%noise_var=input('Enter  noise  variance:  ') 
AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
b=input('Enter  number  of  bits:  ') 
delta=2A-b; 

for  outer_loop=  1:1000 
outeMoop 

f=zeros(2,N); 
xterml=zeros(l,N); 
xterm2=zeros(  1  ,N) ; 
gterml  =zeros(  1  ,N) ; 
gterm2=zeros(  1  ,N) ; 


%  Construct  test  signal  and  interpolation  functions  sampled  at  1  Hz 

for  n=l:N 

arg=(n-N/2)*T; 

f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2))  +  delta*(rand(l,l)-.5); 
f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2))  +  delta*(rand(l,l)-.5); 

xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 

xterm2(n)=(A0-Al*exp(j*2*pi*W*(arg-tau/2)))*f(2,n)/(A0A2-AlA2); 

gterml  (n)=W*sin(pi*W*(-arg-tau/2))*exp(j*pi*W*(-arg-tau/2))/(pi*W*(-arg-tau/2)); 

gterm2(n)=W*sin(pi*W*(-arg+tau/2))*exp(j*pi*W*(-arg+tau/2))/(pi*W*(-arg+tau/2)); 

end 

%for  n=(N/2)+l:N 


%  arg=(n-249)*T; 


%  f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2)); 

%  f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2)); 

%  xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 

%  xterm2(n)=(A0- A 1  *exp(j  *2  *pi*  W*  (arg-tau/2)))  *f(2,n)/(A0A2- A 1 A2) ; 

%  gterml(n)=W*sin(pi*W*(-arg+tau/2))*exp(j*pi*W*(-arg+tau/2))/(pi*W*(- 
arg+tau/2)); 

%  gterm2(n)=W*sin(pi*W*(-arg-tau/2))*exp(j*pi*W*(-arg-tau/2))/(pi*W*(-arg-tau/2)); 
%end 


%f(N/2)  =  2*W; 

%xtennl  (N/2)=(  AO- A 1 )  *f  (N/2)/(A0A2- A 1 A2) ; 
%xterm2(N/2)=(A0-Al)*f(N/2)/(A0A2-AlA2); 
%gterml(N/2)=W; 

%gterm2(N/2)=W; 


%  Interpolate  at .  1  Hz  sampling 

f_hat_0=0; 

for  m=l:N 

f_hat_0  =  f_hat_0  +  xterml(m)*gterml(m)  +  xterm2(m)*gterm2(m); 
end 

f_hat_0_save(outer_loop)=2*real(f_hat_0); 

2*real(f_hat_0); 

%f(l,N/2) 

%f(2,N/2) 

end 


mean(f_hat_0_save) 

var(f_hat_0_save) 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  Sine  function  —  Quantization  Noise  Sims 
%  BW  =  .1  Hz,  Sampled  at  .1  Hz 


Fs=.l 

N=500 

Fc=.  1 

W=.l 

T=l/W 

tau=.l 

%noise_var=input('Enter  noise  variance:  ') 
AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
b=input('Enter  number  of  bits:  ') 
delta=2A-b; 

for  outer_loop=  1:1000 
outeMoop 

f=zeros(2,N); 
xterml=zeros(l,N); 
xterm2=zeros(  1  ,N) ; 
gterml  =zeros(  1  ,N) ; 
gterm2=zeros(  1  ,N) ; 


%  Construct  test  signal  and  interpolation  functions  sampled  at  1  Hz 

for  n=l:N 

arg=(n-N/2)*T; 

f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2))  +  delta*(rand(l,l)-.5); 
f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2))  +  delta*(rand(l,l)-.5); 

xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 

xterm2(n)=(A0-Al*exp(j*2*pi*W*(arg-tau/2)))*f(2,n)/(A0A2-AlA2); 

gterml  (n)=W*sin(pi*W*(-arg-tau/2))*exp(j*pi*W*(-arg-tau/2))/(pi*W*(-arg-tau/2)); 

gterm2(n)=W*sin(pi*W*(-arg+tau/2))*exp(j*pi*W*(-arg+tau/2))/(pi*W*(-arg+tau/2)); 

end 

%for  n=(N/2)+l:N 


%  arg=(n-249)*T; 


%  f(l,n)  =  2*W*sin(2*pi*W*(arg+tau/2))/(2*pi*W*(arg+tau/2)); 

%  f(2,n)  =  2*W*sin(2*pi*W*(arg-tau/2))/(2*pi*W*(arg-tau/2)); 

%  xterml(n)=(A0-Al*exp(j*2*pi*W*(arg+tau/2)))*f(l,n)/(A0A2-AlA2); 

%  xterm2(n)=(A0- A 1  *exp(j  *2  *pi*  W*  (arg-tau/2)))  *f(2,n)/(A0A2- A 1 A2) ; 

%  gterml(n)=W*sin(pi*W*(-arg+tau/2))*exp(j*pi*W*(-arg+tau/2))/(pi*W*(- 
arg+tau/2)); 

%  gterm2(n)=W*sin(pi*W*(-arg-tau/2))*exp(j*pi*W*(-arg-tau/2))/(pi*W*(-arg-tau/2)); 
%end 


%f(N/2)  =  2*W; 

%xtennl  (N/2)=(  AO- A 1 )  *f  (N/2)/(A0A2- A 1 A2) ; 
%xterm2(N/2)=(A0-Al)*f(N/2)/(A0A2-AlA2); 
%gterml(N/2)=W; 

%gterm2(N/2)=W; 


%  Interpolate  at .  1  Hz  sampling 

f_hat_0=0; 

for  m=l:N 

f_hat_0  =  f_hat_0  +  xterml(m)*gterml(m)  +  xterm2(m)*gterm2(m); 
end 

f_hat_0_save(outer_loop)=2*real(f_hat_0); 

2*real(f_hat_0); 

%f(l,N/2) 

%f(2,N/2) 

end 


mean(f_hat_0_save) 

var(f_hat_0_save) 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  Low  Pass  Noise  function 
%  BW  =  .1  Hz,  Sampled  at  .1  Hz 

elf 

clear 


Fs=l 

N=4097 

W=input( 'Enter  W:  ') 

T=fix(l/W) 

%tau=input('Enter  tau:  ') 
tau=2 

%noise_var=input('Enter  noise  variance:  ') 

noise_var=0 

AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
t=-5:5; 

whitenoise=randn(  1 , 1 0000) ; 

F=[0  2*.06  2*.l  1]; 

A=[l  10  0]; 
b=fir2(2048 ,F,  A) ; 
lpnoise=filter(b,  1 .0,whitenoise); 

f=zeros(l,N); 
xterml=zeros(l,N); 
xterm2=zeros(  1  ,N) ; 
gterml=zeros(l,N); 
gterm2=zeros(l,N); 

w=window(@hamming,N)'; 

f=lpnoise(7000-2048:7000+2048); 

%  Construct  test  signal  and  interpolation  functions  sampled  at  1  Hz 
for  n=l:Fs:2048 

xterml(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

xterm2(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

gterml(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

gterm2(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

end 


for  n=2050:Fs:4097 


xterml(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

xterm2(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

gterml(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

gterm2(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

end 


xterm 1  ( 2049)=(  AO- A 1 )  *f(2049)/(A0A2- A 1 A2) ; 
xterm2(2049)=(  AO- A 1 )  *f(2049)/(A0A2- A 1 A2) ; 
gterml(2049)=W; 
gterm2 ( 2049)= W; 


%Interpolate  at  .1  Hz  sampling 

f_hat_save=zeros(l,length(t)); 

for  t_index=l :  length  (t) 

f_hat=zeros(l  ,length(t)); 

%for  m  =  T+length(t)+l:T:N-T-length(t) 
count=0; 

for  m  =  2049-2000:T:2049+2000 

count=count+l; 

m_save(count)=m; 

f_hat(t_index)  =  f_hat(t_index)  +  xterml(m+tau/2)*gterml(t(t_index)+m+tau/2) 
xterm2(  m-tau/2)  *  gterm2  (t(t_index)+m-tau/2) ; 

end 

f_hat_save(t_index)=2*real(f_hat(t_index)); 

end 

f_test=f(2039-fix(length(t)/2):2039+fix(length(t)/2)) 
f  hat  save 


plot(t,f_test) 
hold  on 

plot(t,f_hat_save,'.') 

grid 

xlabel('Time  (Seconds)') 
ylabel('f(t),  fhat(t)') 

title('16-QAM  Eg,  W=.2  Hz,  Fs=.2  Hz,  (Solid  —  Original,  Dots  —  Reconstructed)') 


%  Reconstruction  of  1/2  Nyquist  Rate  Sampled  Signal 
%  Test  Signal  is  the  .25  SQRTCOS  filtered  16-QAM 
%  BW  =  .2  Hz,  Sampled  at  .2  Hz 

elf 

clear 


Fs=l 

N=4097 

W=input( 'Enter  W:  ') 

T=fix(l/W) 

%tau=input('Enter  tau:  ') 
tau=2 

noise_var=input('Enter  noise  variance:  ') 
white_noise=sqrt(noise_var)  *randn(  1  ,N) ; 

AO  =  2/T 

A1  =  (2/T)*cos(pi*W*tau) 
t=-50:50; 


load  -MAT  snr_inf.bak 
qam_sig=zeros(  1  ,N) ; 

qam_sig(2049-1250:2049+1250)=real(MatrixData(100:2600)); 

w=hamming(l,2501); 

%qam_sig(2049-1250:2049+1250)=w.*qam_sig(2049- 1250:2049+1250); 

f=zeros(l,N); 

xterml=zeros(l,N); 

xterm2=zeros(l,N); 

gterml=zeros(l,N); 

gterm2=zeros(l,N); 

f  =  qam_sig  +  white_noise; 

%  Construct  interpolation  functions  sampled  at  1  Hz 
for  n=l:Fs:2048 

xterml(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

xterm2(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

gterml(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

gterm2(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 


end 


for  n=2050:Fs:4097 


xterml(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

xterm2(n)=(A0-Al*exp(j*2*pi*W*(n-2049)))*f(n)/(A0A2-AlA2); 

gterml(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

gterm2(n)=W*sin(pi*W*(2049-n))*exp(j*pi*W*(2049-n))/(pi*W*(2049-n)); 

end 

xterm 1  ( 2049)=(  AO- A 1 )  *f(2049)/(A0A2- A 1 A2) ; 
xterm2(2049)=(  AO- A 1 )  *f(2049)/(A0A2- A 1 A2) ; 
gterml(2049)=W; 
gterm2(2049)=W; 

%Interpolate  at  .2  Hz  sampling 
f_hat_save=zeros(l,length(t)); 
for  t_index=l :  length  (t) 
f_hat=zeros(l  ,length(t)); 

%for  m  =  T+length(t)+l:T:N-T-length(t) 
count=0; 

for  m  =  2049- 1250:T:2049+ 1250 

count=count+l; 

m_save(count)=m; 

f_hat(t_index)  =  f_hat(t_index)  +  xterml(m+tau/2)*gterml(-t(t_index)+m+tau/2)  + 
xterm2  ( m-tau/2)  *  gterm2  (-t(t_index)+m-tau/2) ; 

end 

f_hat_save(t_index)=2*real(f_hat(t_index)); 

end 


f_test=f(2049-fix(length(t)/2):2049+fix(length(t)/2)) 


f_hat_save 

snr=10*logl0(sum(qam_sig.A2)/noise_var) 

mse=sum((f_test-f_hat_save).A2) 

y=2*W*real(conv(xterml,gterml)  +  conv(xterm2,gterm2)); 
y(4097-fix(length(t)/2):4097+fix(length(t)/2)) 

plot(t,f_test) 
hold  on 

plot(t,fJiat_save,'.') 

grid 

xlabel('Time  (Seconds)') 
ylabel('f(t),  fhat(t)') 

title('16-QAM  Eg,  W=.2  Hz,  Fs=.2  Hz,  (Solid  —  Original,  Dots  —  Reconstructed)') 


Appendix  D:  Random  Process  Reconstruction  Convergence  Proof  in  a  Mean 
Squared  Error  Sense 

Let  fit )  be  a  zero  mean  stationary  random  process,  and  fit)  be  the  reconstructed 
random  process  using  the  algorithms  derived  in  the  technical  report  found  in  Appendix  A. 

Although  numerous  simulations  indicate  that  the  reconstruction  equations  are  valid  for  a 
random  process  we  also  desire  to  show  the  reconstruction  error  is  zero  in  a  mean  square 
sense.  That  is  we  desire 


or 


E{f2(t)}-2E{f(t)f(t)}+E{f2(t)}=0. 

Indeed  we  shall  find  that 

E{f\t)}=  E{f(t)f(t)}=  E{f2(t)}=  Rf  (0) 

so  that  the  reconstruction  mean  square  error  is  equal  to  zero.  We  consider  the  three 
expected  value  functions  in  order. 

1.)  First  we  consider  £’{/2(t)}- 

By  definition 

Rf(T)  =  E{f(t)f(t  +  T)} 
and  with  r  =  0  we  find  trivially  that 

Rf(0)  =  E{f2(t)}. 


2.)  Secondly  we  consider  E{f(t)  fit)]. 

It  is  useful  to  recall  that  an  autocorrelation  function  can  be  represented  as  a  sampling 
expansion  since  we  assume  the  autocorrelation  function  is  bandlimited  and  the  random 
process  is  stationary. 

Thus  we  may  express  the  autocorrelation  in  sampled  form  as 


or  setting  t  =  0  we  have 


n=+ oo 

Rf(r)=  ^Rf(nT) 


n=— oo 


sin  \l7lW  {r  -  nT)\ 
2  7tW(T-nT) 


n=+ oo 

*/(°)  =  HRf(nT) 


n=— oo 


sin  (-nT)] 
2  7tW(-nT) 


We  may  write  the  reconstruction  equation  for  /'(f)  as 


/(f)  =  2 Re  £ 


A  ,'2™+r/2)  Wfsin[^_nr  +  T/2)]  7W((_,ir+r/2) 

— - - - - tynl+T/Z) - eJ 

AS -A?  7iW (t  —  nT  +  f/2) 


+ 


A 


0 


Agj2*,,,r"f<2)  f (nT  tI2)W sin -nT-r/2)\ ^ww/2) 
A02  -  A2  7  mit-nT  -t/2) 


Note  that 


£{/(0/(»^  +  /- T  / 2)}  =  //  (f  -nT-l+tl 2)  =  (-nr  -l+r/2) 

since,  assuming  stationarity,  we  may  set  t  =  0  for  convenience. 

Since  /(f)  is  not  a  function  of  summation  index  n  we  may  move  the  term  inside  the 
summation  when,  using  the  preceding  equation  and  after  moving  the  expectation  operator 
inside  the  summation  as  well  (expectation  is  a  linear  operation),  we  form 


r  v  n=+ 00 

£l/W(0/=2Re  X  Rf(-nT-T/ 2) 


An  -  A,e 


j2m(nT+zl2) 


Wsin[nW(-nT  +  T/2)] 


j7tW(-nT+Tl2 ) 


A2 -A2 


7tW{—nT  +  t/2) 


+ 


R  r  (—nT  +  T  /  2)  - 


-  Aiejm-T,2)  W sin [7tW{-nT  -t/2)\ 
A]  -  A2  7iW(-nT-rl 2) 


Note  that  this  expression  reconstructs  Rf  (O)from  a  Vi  rate  sampled  version  of  the 
autocorrelation  function,  thus  we  have 

E{f(t)f(t)}=Rf( 0). 


3.)  We  abbreviate  the  final  exercise  showing  E\f2  (t)}=  Rf  (0)  since  the  development 

parallels  the  preceding  derivation  which  found  that  £’{/(t)/(t)}=  R ,  (0) .  Additionally, 

the  manipulations  are  very  tedious  due  to  the  product  of  the  two  reconstruction 
summations.  The  key  steps  are  to:  (i)  Replace  the  product  of  sums  by  a  double 
summation  of  products  using  two  summation  indices  such  as  n  and  k  to  keep  the  terms 
separate,  (ii)  expand  the  four  terms  created,  (iii)  move  the  expectation  operator  inside  the 
double  summation  to  form  the  autocorrelation  functions  of  the  form  Rf  (n  -k),  (iv)  then 

finally  observe  that  the  sinusoid  products  are  orthogonal  and  sum  to  zero  over  infinity 
except  for  n=k  leading  to  R ,  (0) . 


Then  we  find 


E{f2(t)}=Rf  (0). 


